%%%% equilibrium outcomes with constrained solution

function [K_1_grid, tau_k_2_grid, L_1_grid, L_2_grid, q_1_grid, B_0_grid, tau_l_2_grid] = eqm_constrained(B_e_0, par, grid_number, PSI_step)

K_1_grid     = zeros(grid_number,1); 
L_1_grid     = zeros(grid_number,1); 
L_2_grid     = zeros(grid_number,1); 
tau_k_2_grid   = zeros(grid_number,1);
tau_l_2_grid   = zeros(grid_number,1);
q_1_grid     = zeros(grid_number,1);
B_0_grid     = zeros(grid_number,1);


for jj = 1:1:grid_number
    PSI     = PSI_step * (jj - 1);
    
    K_1 = ((1 + PSI) / (1 + (1 + par.NU) * PSI) * par.A * (1 - par.ALPHA) / par.MU)^(1 / par.NU) *...
    ( ((par.phi + PSI) / par.phi / (1 + PSI) / par.BETA - (1 - par.DELTA)) / par.A / par.ALPHA )^((par.ALPHA + par.NU) / par.NU / (par.ALPHA - 1));

    K_1_L_2 = ((par.phi + PSI) / par.phi / (1 + PSI) / par.BETA / par.ALPHA)^(1 / (par.ALPHA - 1));
    
    L_1 = (1 + PSI) / (1 + (1 + par.NU) * PSI) * par.A  / par.MU;
    %L_2 = (1 + PSI) / (1 + (1 + par.NU) * PSI) * par.A * (1 - par.ALPHA) * K_1_L_2^par.ALPHA / par.MU;
    L_2 = (par.A * (1 - par.ALPHA) / par.MU * (1 + PSI) / (1 + (1 + par.NU) * PSI))^(1 / (par.ALPHA + par.NU))...
        * K_1^(par.ALPHA / (par.ALPHA + par.NU));
    q_1     = (K_1 - B_e_0) / par.phi / K_1;
    r_2     = par.A * par.ALPHA * (K_1 / L_2)^(par.ALPHA - 1) + 1 - par.DELTA;
    tau_k_2   = 1 - q_1 / par.BETA / r_2;
    w = par.A * (1 - par.ALPHA) * (K_1 / L_2)^par.ALPHA;
    tau_l_2 = 1 - L_2^par.NU * par.MU / w;

    C_1     = par.A * L_1 - K_1 - par.G_1; % goods mkt clearing in period 1
    C_2     = par.A * K_1^par.ALPHA * L_2^(1 - par.ALPHA) - par.G_2;
    B_0     = (C_1 - par.MU * L_1^par.NU * L_1) + par.BETA * (C_2 - par.MU * L_2^par.NU * L_2) - (1 / par.phi - 1) * max(0, K_1 - B_e_0 / (1 - par.phi));

    K_1_grid(jj) = K_1;
    L_1_grid(jj) = L_1;
    L_2_grid(jj) = L_2;
    tau_k_2_grid(jj) = tau_k_2;
    tau_l_2_grid(jj) = tau_l_2;
    q_1_grid(jj) = q_1;
    B_0_grid(jj) = B_0;
end

end